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ABSTRACT 

We present a comprehensive investigation of the cosmological evolution of the luminosity 
function of galaxies and active galactic nuclei ( AGN) in the infrared (IR) . Based on the observed 
dichotomy in the ages of stellar populations of early-type galaxies on one side and late- type galax- 
ies on the other, the models interprets the epoch-dependent luminosity functions at z > 1.5 using 
a physical model for the evolution of proto-spheroidal galaxies and of the associated AGNs, while 
IR galaxies at z < 1.5 are interpreted as being mostly late- type "cold" (normal) and "warm" 
(starburst) galaxies. As for proto-spheroids, in addition to the epoch-dependent luminosity func- 
tions of stellar and AGN components separately, we have worked out, for the first time, the 
evolving luminosity functions of these objects as a whole (stellar plus AGN component), taking 
into account in a self-consistent way the variation with galactic age of the global SED. The model 
provides a physical explanation for the observed positive evolution of both galaxies and AGNs 
up to z ~ 2.5 and for the negative evolution at higher redshifts, for the sharp transition from 
Euclidean to extremely steep counts at (sub-)mm wavelengths, as well as the (sub-)mm counts 
of strongly lensed galaxies, that are hard to account for by alternative, physical or phenomeno- 
logical, approaches. The evolution of late-type galaxies and of z < 1.5 AGNs is described using 
a parametric phenomenological approach. The modeled AGN contributions to the counts and 
to the cosmic infrared background (CIB) are always subdominant. They are maximal at mid-IR 
wavelengths: the contribution to the 15 and 24 /im counts reaches 20% above 10 and 2 mJy, 
respectively, while the contributions to the CIB are of 8.6% and of 8.1% at 15 /xm and 24 fim, 
respectively. The model provides a good fit to the multi-wavelength (from the mid-IR to mil- 
limeter waves) data on luminosity functions at different redshifts and on number counts (both 
global and per redshift slices). A prediction of the present model, useful to test it, is a systematic 
variation with wavelength of the populations dominating the counts and the contributions to the 
CIB intensity. This implies a specific trend for cross- wavelength CIB power spectra, that is found 
to be in good agreement with the data. 
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Subject headings: galaxies: formation - galaxies: evolution - galaxies: elliptical - galaxies: high 
rcdshift - submillimeter 



Introduction 



The huge amount of infrared (IR) to millimeter- wave data that has been accumulating in the last several 
years have not yet led to a fully coherent, established picture of the cosmic star-formation history, of the IR 
evolution of Active Galactic Nuclei (AGNs), and of the inter-relations between star formation and nuclear 
activity. 



and AGN luminosity functions over a broad wavelength ranee have been worked out (e.g. Bcthermin et al. 


2012a 


. 2011: GruDDioni et al. Iboil 


; Rahmati & van der Werf 12011; Marsden et al. Iboil 


Franceschini et al. 


2010; 


Valiante et al. 2009; Le Borgne et al. 20091 Rowan- Robinson 2009). These models generally include 



multiple galaxy populations, with different spectral energy distributions (SEDs) and different evolutionary 
properties, described by simple analytic formulae. In some cases also AGNs are taken into account. All of 
them, however, admittedly have limitations. 

The complex combination of source properties (both in terms of the mixture of SEDs and of evolu- 
tionary properties), called for by the richness of data, results in a large number of parameters, implying 
substantial degeneracies that hamper the interpretation of the results. The lack of constraints coming from 
the understanding of the astrophysical processes controlling the evolution and the SEDs limits the predictive 
capabilities of these models. In fact, predictions of prc-Herschel phenomenological models, matching the 
data then available, yielded predictions for Herschel counts quite discrepant from each other and with the 
data. 

The final goal is a physical model linking the galaxy and AGN formation and evolution to primordial 
density perturbations. In this paper we make a step in this direction presenting a comprehensive 'hybrid' 
approach, combining a physical, forward model for spheroidal galaxies and the early evolution of the associ- 
ated AGNs with a phenomenological backward model for late-type galaxies and for the later AGN evolution. 
We start from the consideration of the observed dichotomy in the ages of stellar populations of early-type 
galaxies on one side and late- type galaxies on the other. Early-type galaxies and massive bulges of Sa galaxies 
are composed of relatively old stellar populations with mass-weighted ages of > 8-9 Gyr (corresponding to 
formation redshifts z > 1—1.5), while the disc components of spiral and irregular galaxies are characterized 
by significantly younger stellar populations. For instance, the luminosity-weighted age for most of Sb or 
later- type spirals is < 7 Gyr (cf. Bernardi et al. 2010, their Fig. 10), corresponding to a formation redshift 
z < 1. Thus proto-spheroidal galaxies are the dominant star-forming population at z > 1.5, while IR galaxies 
at z < 1.5 are mostly late-type "cold" (normal) and "warm" (starburst) galaxies. 

Fuller hierarchical galaxy formation models, whereby the mass assembly of galaxies is related to structure 
formation in the dark matter and the star formation and merger histories of galaxies of all morpho logical types 



are calculated based on physical prescriptions have been recently presen ted by several groups (jLacey et al 



2008 ; Fontanot et al. 2009 ; Narayanan et al. 20ld Shimizu et al. 2012 ). However, the predictions for the 
IR evolution of galaxies are limited to a small set of wavelengths and frequently highlight serious diff iculties 
with accounting for observational data ( Lacev et al. Il2010l : iNiemi et al. 11201 2l : lHavward et al. 112012 ). 
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While the evolution of dark matter halos in the framework of the 'concordance' ACDM cosmology is rea- 
sonably well underst ood thanks to N-body simulations such as the Millennium, the Millennium- XXL and the 
Bolshoi simulations ( Springel et al. 20051 Bovlan-Kolchin et al. 2009 ; Angulo et al. 2012 ; Klypin et al 



20111 ). establishing a clear connection between dark matter halos and visible objects proved to be quite 
challenging, especially at (sub-)mm wavelengths. The early predictions of the currently favoured scenario, 
whereby both the star-formation and the nuclear activity a re driven by mergers , were more than one order 

The 



of magnitude below the observed SCUBA 850 fim counts (jKaviani et al 



2003 



Baugh et al. 



200 



5). 



basic problem is that the duration of the star- formation activity triggered by mergers is too short, requir- 
ing non standard assumptions either on the Initial Mass Function (IMF) or on dust properties to account 
for the measured source counts. The problem is more clearly illustrate d in terms of redshift-dependent far 
IR/sub-mm luminosity function, estimated on the basis of Herschel data ( Eales et al. Iboiol Gruppioni et al. 



20ld lLapi et al 
SFR ~ 300 M Q yr 



20111 ). These estimates consistently show that z ~ 2 galaxies with Star Formation Rates 
_1 have comoving densities $300 ~ 10~ 4 Mpc~ 3 dex -1 . The comoving density of the cor- 
responding halos is n(A:f v i r ) ~ ^soo^oxp/tsfr), where M v ; r is the total virial mass (mostly dark matter), 
tsfr is the lifetime of the sta r-forming phase an d t exp is the expansion timescale. For the fiducial lifetime 

12. 92 while for tsfr p 0.1 Gyr, typical 
(|2011l ) model implies a 



tsfr — 0.7 Gyr advocated by lLapi et al 



j201lD, log(M vir /M ) 
i/Mq) ~ 12.12. Thus while the lLapi et al 



of a merger-driven starburst, log(M v 
SFR/Af v ; r ratio easily accounted for on the basis of standard IMFs and dust properties, the latter scenario 
requires a SFR/Af v ; r ratio more than a factor of 6 higher. 



To reach the required values of SFR/Af v j r or , equivalently, of Lnt/Mvir, iBaugh et al. ( 2005 ) resorted 



to a top-heavy IMF while iKaviani et al. I (|2003l ) assumed that the bulk of the sub-mm emission comes 



from a huge amount of cool dust. But even tweaking with t he IMF and with dust properties, fit s of 



the sub-mm counts obtained within the merger-driven scenario ([Lacev et al. luOlOt iNiemi et al. 112012 ) are 
generally unsatisfactory. Further constraints on physical models come from the clustering propert ies of 



sub-mm galaxies that are determined by their effective halo masses. As shown by Xia et al. ( 2012h both 
the angular correlation function of detected sub-mm galaxies and the power spectrum of fluctuations of the 
cosmic infrared background indicate halo masses larger than implied by the major mergers plus top-heavy 



initial stellar mass function scenario (|Kim et al. 1 1201 lh and smaller than imp l ied by cold flow models but 



consistent with the self-regulated baryon collapse scenario (|Granato et al. 112004 lLapi et al. 200G; lLapi et al. 
201ll ). 



As is well known, the strongly negative K-correction emphasizes high-z sources at (sub-)mm wave- 
lengths. The data show that the steeply rising portion of the (sub-)mm counts i s indeed dominated by 



ultra-luminous star-fo r ming galaxies with a redshift distributi on peaking at z c± 2.5 (Chapman et al. 20051 : 



Aretxaga et al. 1 120071 : lYun et al. 1 12012b ISmolcic et al. 1120121 ). As shown bv lLapi et al. I (|2011l ). the self- 



rcgulated baryon collapse scenario provides a good fit of the (sub-) mm data (counts, redshift-dependent 
luminosity functions) as well as of the stellar mass functions at dif ferent redshifts. Moreover, the counts of 
strongly lensed galaxies were pred icted with remarkable accuracy jNegrello et al. 1120071 lioict lLapi et al 



2012t iGonzalez-Nuevo et al. II2012T). Furthe r considering that this scenario accounts for the clustering prop- 



erties of sub-mm galaxies (|Xia et al 

present analysis. However, we upgrade this model in two respects 



20121 ) . we conclude that it is well grounded, and we adopt it for the 

First, while, on one side, the model 



envisages a co-evolution of spheroidal galaxies and active nuclei at their centers, the emissions of the two 
components have been, so far, treated independently of each other. This is not a problem in the wavelength 
ranges where one of the two components dominates, as in the (sub-)mm region where the emission is dom- 
inated by star-formation, but is no longer adequate at mid-IR wavelengths, where the AGN contribution 
may be substantial. In this paper we present and exploit a consistent treatment of proto-spheroidal galaxies 
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including both components. Second, while the steeply rising portion of (sub-)mm counts is fully accounted 
for by proto-spheroidal galaxies, late-type (normal and starburst) galaxies dominate both at brighter and 
fainter flux densities and over broad flux density ranges at mid-IR wavelengths. At these wavelengths, AGNs 
not associated to proto-spheroidal galaxies but either to evolved early type galaxies or to late-type galax- 
ies are also important. Since we do not have a physical evolutionary model for late- type galaxies and the 
associated AGNs, these source populations have been dealt with adopting a phenomenological approach. 

Another distinctive feature of the present model is that we have attempted to fit simultaneously the 
data over a broad wavelength range, from mid-IR to mm waves. As mentioned in several papers, this faces 
us with several challenges. First, the data co me from different instru ments and the relative calibration is 



sometimes problematic (see the discussion in iBethermin et al. _ 2011 ) . Systemat ic calibration offsets may 



hinder simultaneous fits of different data sets. For example, iMarsden et al. I (|201lh pointed out that there is 
considerable tension between the SCUBA 850 /im counts and the AzTEC counts at 1.1mm, and indeed the 
850 fim and mm-wave counts have been repeatedly corrected (generally downwards) as biases were discovered 
and better data were acquired. Also, the very complex SEDs in the mid-IR, where strong polycyclic aromatic 
hydrocarbon (PAH) emission features show up, make the counts exceedingly sensitive to the details of the 
spectral response function of the specific instrument and introduce large uncertainties in the conversion from 
broad-band measurements to monochromatic flux densities giving rise to strong discrepancies among data 
sets nominally referring to the same wavelength. In fact, large discrepancies are present among different 
determinations of 15 /im and 60 /im source counts. 

The plan of the work is the following. In Section[2]we describe the physical model for the evolution of 
proto-spheroidal galaxies and of the associated AGNs and the SEDs adopted for these sources. Section |3] 
deals with the evolutionary model for late- type galaxies and z < 1.5 AGNs. In Section[3]we present the 
formalism to compute the source counts of unlensed and lensed sources, the cumulative flux density as a 
function of redshift and the contributions to the CIB. In Section[5]we report on the determination of the 
best fit values of the model parameters. In Section[6] the model results are compared with data on multi- 
frequency luminosity functions at various redshifts and on source counts, both total and per redshift slices. 
The multi-frequency power-spectra of CIB fluctuations implied by the model are discussed in Section[7] 
Finally, Section[8] contains a summary of the paper and our main conclusions. 

Tabulations of multi-frequency model counts, redshift distributions, SEDs, redshift-dependent luminos- 
ity functions at several wavelengths, and a large set of figures comparing model predictions with the data 



are available in the Web site http://people.sissa.it/~zcai/galaxy_agn/ 



Throughout this paper we adopt a flat cosmology with present day matter and baryon density, in units 
of the critical density, r2 m = 0.27 and f?b o = 0.044; Hubble constant h = H o /100 = 0.71; spectrum of 
primordial density perturbations with slope n = 1 and normalization on a scale of 8h _1 Mpc a s = 0.81. 



2. Star-forming proto-spheroidal galaxies 
2.1. Overview of the model 



We adopt the model by iGranato et al 



(|2004 see also Lapi et al. 2006, 2011; Mao et al. 2007) that 
interprets powerful high-z sub-mm galaxies as massive proto-spheroidal galaxies in the process of forming 
most of their stellar mass. It hinges upon high resolution numerical simulations showing t hat dark matter 
halos form in two stages (jZhao et al. I l2003t IWang et al. 1 1201 lc lLapi fc Cavaliere 1 1201 If ) . An early fast 
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collapse of the halo bulk, including a few major merger events, reshuffles the gravitational potential and 
causes the dark matter and stellar components to undergo (incomplete) dynamical relaxation. A slow 
growth of the halo outskirts in the form of many minor mergers and diffuse accretion follows; this second 
stage has little effect on the inner potential well where the visible galaxy resides. 

The star formation is triggered by the fast collapse/merger phase of the halo and is controlled by self- 
regulated baryonic processes. It is driven by the rapid cooling of the gas within a region with radius ~ 30% 
of the halo virial radius, i.e. of ~ 70(M v i r /10 13 M Q ) 1 / 3 [(1 + z v j r )/3]~ 1 kpc, where M v - lr is the halo mass and 
z v i r is the virialization redshift, encompassing about 40% of the total mass (dark matter plus baryons). The 
star formation and the growth of the central black- hole are regulated by the energy feedback from supernovae 
(SNe) and from the active nucleus, is very soon obscured by dust and is stopped by quasar feedback. The 
AGN feedback is relevant especially in the most massive galaxies and is responsible for their shorter duration 
(5 — 7 x I0 8 yr) of the active star-forming phase. In less massive proto-spheroidal galaxies the star formation 
rate is mostly regulated by SN feedback and continues for a few Gyr. Only a minor fraction of the gas 
initially associated to the dark matter halo is converted into stars. The rest is ejected by feedback processes. 

The equations governing the evolution of the baryonic matter in dark matter halos and the adopted 
values for the parameters are given in the Appendix where some examples of the evolution with galactic 
age (from the virialization time) of quantities related to the stellar and to the AGN component are also 
s hown. For additiona l details and estim ates of physically plaus ible ranges for each parameter we refer to 



Granato et al. I (|2004) . lLapi et alJ (|2006h and 



be in passive evolution at z < 1 
higher redshifts. 



Mao et al 



(|2007fl . Since spheroidal galaxies are observed to 



1.5 (e.g. iRenzini II2006I ). they are bright at sub-mm wavelengths only at 



2.2. Luminosity functions 



The bolometric luminosity function (LF) of proto-spheroids is obtained convolving the halo formation 
rate dNsT(M v [ r , z)/dt with the galaxy luminosity distribution, P(L, z; M v - U ). The halo formation rate is well 
approximated, for z > 1.5, by the positive term of the cosmic time derivative of the halo mass function 
iVs T- For the latter, giving the average comoving number density of haloes of given mass, M v ; r , we adopt 
the lSheth fc Tormen I |l999h analytical expression 



N S T{M vir ,z)dM viT 



Pm,0 t i 



d\vtv 
dlni\fvi 



■dM viT 



(1) 



where p m .o = fl m ,oPc,o is the present day mean comoving matter density of the universe and v = [5 c (z) / a(M v i I )] 2 , 
with 5 c (z) = 8o(z)D(0)/ D(z). The critical value of the initial overdensity t hat is required for spherical col- 
lapse at z, S (z), is 8 c (z) = 8 Q (z)D(0)/D(z) with (jNakamura fc Suto Ifl997f) 



So(z) 



3(12tt) 



2/3 



20 



-[1 + 0.0123 log Q m (z)] ~ 1.6865[l + 0.0123 log Sl m (z)]. 



1991 



Carroll et al. 



199 



2) 



The linear growth factor can be approximated as (|Lahav et al. 

w 2(1 + z)I [70 140 y ' 140 mV ' m y> 

The mass variance a(M v - lr ) of the primordial perturbation fi eld smoothed o n a sc ale containing a mass 
M V - 1T with a top-hat window function was computed using the iBardeen et al. I (|1986T ) power spectrum with 
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correction for baryons (jSugivama I Il995). for our choice of cosmological parameters (see Section[TJ. The 
results are accurately approximated (error < 1% over a broad range of M v j r , 10 6 < M v ; r /M < 10 16 ) by 



a(M vir ) 



0.8 
(184 



[14.110393 - 1.1605397a; - 0.0022104939x 2 



+ 0.0013317476a; 3 - 2.1049631 x lO^a; 4 ] 
where x = log(M v i r /M ). Furthermore 

f ST (v)=A[l + (av)-v}( Y ) — rjr , 

where A = 0.322, p = 0.3, and a = 0.707. 
The halo formation rate is then 

dN ST (M vil ,z) 



(2) 



dt 

-N ST {M vil ,z) 



7V ST (M vir ,z) 



rfln/ ST (t/) 
dt 



aS c 2p a 



■2p 



(J 2 S c a ^v + a P5l p S c 



dS c dz 
dz dt 



N ST (M viT ,z) 



P 



av 

T + 1 + {av)P 



c?ln v 



dz 



dz 
It 



(3) 



where dz/dt = -H (l + z)E(z) with E(z) = ^/Oa,o + ^m,o(l + z) 3 . 



The comoving differential luminosity function $(log L, z), i.e. the number density of galaxies per unit 
logL interval at redshift z, is given by 



$(logi,z) 



M max 
M min 

P(logL,z;M 



dz v 



dt v 



dz v 



dN ST 



dt v 



(M, 



viri ^vir 



vm ^vir 



), 



(4) 



where P(logL, z; Af v ; r , z V [ r ) is the luminosity distribution of galaxies at redshift z inside a halo of mass M v i r 
virialized at redshift z w -„. We set z™ n = 1.5 and z™^ x = 12. 

As mentioned in Section ^.ll the total luminosity of a galaxy is the sum of those of the stellar component 
and of the active nucleus. For each component we assume a log-normal luminosity distribution 



F[logL|logL]dlogL = 



exp[- log 2 (£/Z.)/2cr 2 

V27TCT 2 



■dlogL, 



(5) 



with dispersion cr, = 0.10 around the mean stellar luminosity L„(z; M v j r , z V j r ) and <t, = 0.35 around the mean 
AGN luminosity X.(z; M v j r , z v i r ). The mean luminosities are computed solving the equations detailed in the 
Appendix. The higher luminosity dispersion for the AGN component reflects its less direct relationship, 
compared to the stellar compone nt, with M v ; r and z v i r . The distribution of the total luminosity, L to t = 
L» + L., is then (jPufresne Il2004h 

P[log L t ot | log £*, log L.]d log Ltot = dlogL tot 

log -Ltot J T T 

"" expj-^-logZ,) 2 ^ 2 } 



27T(T*CT. Ltot - lO 1 

x exp{-[log(L-tot - 10 x ) 



logL.] 2 ^. 2 }. 



(6) 
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In the upper left panel of Fig. Q] we show, as an example, the bolometric luminosity functions at z = 1.5 of 
the stellar and of the AGN components, as well as the luminosity function of the objects as a whole. As 
shown by eq. ([5]) the latter is different from the sum of the first two, although in this case the difference is 
difficult to perceive. The bright end is dominated by QSOs shining unobstructed after having swept away the 
interstellar medium of the host galaxy. In this phase the QSOs reach their maximum luminosity. Around 
log(Xboi/£o) — 13 the AGNs and the starbursts give similar contributions to the bolometric luminosity 
function. The inflection at log(LiR/i0) ~ 11.7 corresponds to the transition from the regime where the 
feedback is dominated by supernovae (lower halo masses) to the regime where it is dominated by AGNs. 
While the star formation in massive halos is abruptly stopped by the AGN feedback after 0.5-0.7 Gyr, it 
lasts much longer in smaller galaxies, implying a fast increase of their number density. 

The upper right panel of the same figure illustrates the evolution with cosmic time of the global luminos- 
ity function. The cooling and free-fall timescales shorten with increasing redshift because of the increase of 
the matter density and this drives a positive luminosity evolution, thwarted by the decrease in the comoving 
density of massive halos. The two competing factors result, for both the starburst and the AGN component 
(see the lower panels of the figure), in a positive evol ution up to z ~ 2 . 5 follo wed by a decline at hi gher 
z, consistent with the observatio nal determinations by Gruppioni et al. ( 2010l) and Lapi et al. ( 2011 ) for 



the starburst component and by Assef et al. ( 2011 ) and Brown et al. ( 20061) for AGNs. The decrease of 



the luminosity function at low luminosities, more clearly visible at the higher redshifts, is an artifact due to 
the adopted lower limit for the considered halo masses. This part of the luminosity function however does 
not contribute significantly to the observed statistics and therefore is essentially irrelevant here. Below the 
minimum virialization redshift, z™" = 1-5, the bolometric luminosity function of proto-spheroidal galaxies 
rapidly declines as they evolve towards the 'passive' phase. The decline is faster at the bright end (above 
log(Lboi/io) — 12) since the switching off of the star formation for the more massive halos occurs on a 
shorter timescale. 

The monochromatic luminosity functions of each component or of objects as a whole can obviously be 
computed using the same formalism, given the respective spectral energy distributions (SEDs). We define 
C* >v = i>L* iV = vf* (i/JZ^ir, C, >u = vh,, v = vf,{v)L,^ \, and L v = £* <v +£, tV , where f{v) is the normalized 
SED (Jdv f(v) = 1). 

Since the model cannot follow in detail the evolution of the AGN SEDs during the short phase when 
they shine unobstructed by the interstellar medium of the host galaxy, t he distinction betw een obscured 
and unobscured AGNs in the model is made in two ways. First, following Lapi et al. ( 20061) . we choose a 
fixed optical (B-band) "visibility time" , At V [ S = 5x 10 7 yr, consistent with current estimates of the optically 
bright QSO phase. Alternatively, we set the beginning of the optical bright phase at the moment when the 
gas mass fraction is low enough to yield a low optical depth. We estimate that this corresponds to a gas 
fraction within the dark matter potential well / gas = M gas /M V i r < / gas , C rit = 0.03. The two approaches give 
very similar results and we have chosen the criterion / gas < / gas ,crit to compute the luminosity functions at 
optical wavelengths. 



2.3. Spectral energy distributions 



Although there is evidence that the galaxy SEDs vary with luminosity (e.g. lSmith. et al. Il2012l ). lLapi et al. 



(|201ll ) have shown that the sub-mm data can be accurately reproduced using a single SED for proto- 
spheroidal galaxies, i.e. the SED of the strongly lensed z ~ 2.3 galaxy SMM J2135-0102 ([Swinbank et al. 



- 8 - 



20ld: llvison et al. Il201oh . modeled using GRASIL (jsilva et al. 1 1 1 99sh - The basic reason for the higher uni- 



formity of the SEDs of high-z active star-forming galaxies compared to galaxies at low-z is that the far-IR 
emission of the former objects comes almost entirely from dust in molecular clouds, heated by newly formed 
stars, while in low- z galaxies there are important additional contributions from colder 'cirrus' heated by 
older stellar populations. 

This SED worked very well at sub- mm wavelengths but yielded mm- wave counts in excess of the observed 
ones. To overcome this problem the sub-mm slope of the SED has been made somewhat steeper, preserving 
the consistency with the photometric data on SMM J2135-0102 (see Fig. \2§- Moreover, the SED used by 
Lapi et al. has a ratio between the total (8-1000 /im) IR a nd the 8 (im lumino sity (IR8 = LmjL%) of ~ 30, 



far higher than the mean value for z ~ 2 galaxies (IR8 ~ 9. 1Reddv et al. 1120121 ). We have therefore modified 
the near- and mid-IR portions of the SED adopting a shape similar to that of Arp 220. The contribution of 
the passive evolution phase of early-type galaxies is small in the frequency range of interest here and will be 
neglected. 

As mentioned in Section l!01 the model follows the AGN evolution through two phases (a third phase, 
reactivation, will be considered in Section l3~2"|) . For the first phase, when the black- hole growth is enshrouded 
by the abundant, dusty interstellar medium (ISM) of the host galaxy, we ad opt the SED of a heavily 
absorbed AGN taken from the AGN SED library by iGranato fc Danese I (|l994h . Note that these objects 



differ from the classical type-2 AGNs because they are not obscured by a circum-nuclear torus but by the 
more widely distributed dust in the host galaxy. They will be referred to as type-3 AGNs. In the second 
phase t he AGN shines after ha ving swept out the galaxy ISM. For this phase we adopted the mean QSO 



SED bv lRichards et al. I (|2006l ) extended to sub-mm wavelengths assuming a grey-body emission with dust 
temperature Td us t = 80 K and emissivity index /3 = 1.8. These SEDs imply that the IR (8-1000 /Ltm) band 
comprises 92% of the bolometric luminosity of obscured AGNs and 19% of that of the unobscured ones. As 
illustrated by Fig. [21 except in the rare cases in which the AGN bolometric luminosity is much larger than 
that of the starburst, the AGN contribution is small at (sub-)mm wavelengths, while it is important and may 
be dominant, in the mid-IR. This implies that the statistics discussed here are insensitive to the parameters 
describing the extrapolation of the Richards et al. SED to (sub-)mm wavelengths. 

Figure [3] shows the global SEDs and the contributions of the stellar and AGN components for 2 galaxy 
ages and three host halo masses virialized at z v j r = 3. The shorter evolution timescale of the AGNs is clearly 
visible. It is worth noticing that the effect of feedback as a function of halo mass on the SFR is very different 
from that on accretion onto the supermassive black-hole. In the less massive halos the AGN feedback has 
only a moderate effect on the evolution of the SFR and of the accretion rate, that are mostly controlled by 
the SN feedback. With reference to the figure, for \og(M v -„/Mo) — 11.4, the star-formation continues at 
an almost constant rate for a few Gyrs. On the other hand the accretion rate onto the central black-hole 
is at the Eddington limit only up to an age of ~ 0.3 Gyr and afterwards drops to a strongly sub-Eddington 
regime. This is because the growth rate of the reservoir is approximately proportional to the SFR (and 
therefore slowly varying for few Gyrs) while the accretion rate grows exponentially until the mass contained 
in the reservoir is exhausted. From this moment on the accretion rate is essentially equal to the (strongly 
sub-Eddington) inflow rate. For more massive halos the quenching of both the SFR and of the accretion 
occurs more or less simultaneously at ages of ~ 0.5-0.6 Gyr, but while the SFR stops very rapidly, the AGN 
activity continues until the flow of the matter accumulated in the reservoir runs out. At ages £ 0.6 Gyr the 
more massive galaxies are in passive evolution and therefore very weak in the far-IR while star-formation 
and the dust emission are still present in lower-mass galaxies. 
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3. Low redshift (z < 1.5) populations 



3.1. Late-type and starburst galaxies 

We consider two z < 1.5 galaxy populations: "warm" starburst galaxies and "cold" (normal) late- type 
galaxies. For the IR luminosity function of both populations we adopt the functional form advocated by 
Saunders et"aH (jl990h : 



$(logi IR ,z)dlogL m = $* 



L* 



x exp 



\og 2 (l + L m /L*) 



2a 2 



d\og L m 



(J) 



where the characteristic density $* and luminosity L* , the low-luminosity slope a and the dispersion a of each 
population are, in principle, free parameters. However, the low-luminosity portion of the luminosity function 
is dominated by "cold" late-type galaxies and, as a consequence, the value of a of the warm population is 
largely unconstrained; we have fixed it at a warm = 0.01. In turn, the "warm" population dominates at high 
luminosities so that the data only imply an upper limit to u co id- We have set <7 C oid = 0.3. 

For the "warm" population we have assumed power law density and luminosity evolution [<5>*(z) = 
$o(l + z) Q ' 1 '; L*(z) — Lq(1 + z)° l ] up to Zbroak = 1, and »l being free parameters. T he "cold" population 



comprises normal dis c galaxies for which chemo/spectrophotometric evolution models (jMazzei et al. 1 1 1992 : 



Colavitti et al. 20081) indicate a mild (a factor ~ 2 from z = to z = 1) increase in the star formation rate, 
hence of IR luminosity, with look-back time. Based on these results we take, for this population, a.\, = 1 
and no density evolution. At z > Zbrcak both <J>*(z) and L*(z) are kept to the values at Zbrcak multiplied by 
the smooth cut-off function {1 — erf[(z — z cuto ff)/Az)]}/2, with z C utoff = 2 and Az = 0.5. The choice of the 
redshift cutoff for both populations of late-type galaxies is motivated by the fact that the disc component 
of spirals and the irregular galaxies are characterized by relatively young stellar populations (formation 
redshift z < 1—1.5). Above z — 1.5 proto-spheroidal galaxies (including bulges of disk galaxies) dominate the 
contribution to the luminosity function, at least in the observationally constrained luminosity range. The 
other parameters are determined by minimum \ 2 fits to selected data sets, as described in Sect. O Their 
best fit values and the associated uncertainties are listed in Table [TJ 

Although there is cle ar evidence of systematic variations of the IR SEDs of low-z galaxies with luminosity 
(e.g., ISmith et al. 1120121 ). we tried to fit the data with just 2 SEDs, one for the "warm" a nd one for the 
"cold" population. These SEDs were generated by combining those o f iDale fc Helou I (|2002j ). that are best 
determined at mid-IR w avelengths, wit h those of ISmith et al. I (|2012l ). primarily based on Herschel data in 
the range 100-500 ^m. IDale k, Helou I (|2002n give SED templates for several values of the 60 to 100 fj,m 
flux density ratio, log [ f v ( 6 ym ) / f v ( 1 ixm ) | . Usi ng the relation between this ratio and the 3 to 1100 ^m 



luminosity, £tir, given bv lChapman et al. I (|2003l ) we established a correspondence between their SEDs and 



those by Smith et al., labeled by the values of log(LiR,/L0). The combined SEDs are based on Smith et 
al. above 100 /xm and on Dale & Helou at shorter wavelengths. By trial and error we found that the best 
fit to the data is obtained using for the "cold" population the SED corresponding to log(LiR,/L©) = 9.75 
(actually the SEDs change very slightly for log(LrR,/ X Q ) < 9.75) and for the "warm" population the SED 
corresponding to \og(Lm/ Lq) = 11.25. These 2 SEDs are displayed in Fig. 2) 



3.2. Reactivated AGNs 



In the framework of our reference galaxy and AGN evolutionary scenario, most of the growth of super- 
massive black holes is associated to the star forming phase of spheroidal components of galaxies at z > 1.5 
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when the great abundance of interstellar medium favours high accretion rates, at, or even slightly above, 
the Eddington limit. At later cosmic times the nuclei can be reactivated by, e.g., interactions, mergers 
or dynamical instabilities. The accretion rates are generally strongly sub-Eddington. Our evolutionary 
scenario cannot predict their amplitudes and duty cycles. We therefore adopted, also for these objects, a 
phenomenological backward evolution model analogous to that used for the "warm" galaxy population, i.e. 
luminosity functions of the same form of eq. (JTJ) and power-law density and luminosity evolution with the 
same break and cutoff redshifts. However the parameters of the luminosity functions refer to 12 /im (see 
Sect. 13.11 and Tabled]). The data do not allow a determination of the slopes, a, of the faint portions of the 
luminosity functions. We have set a = 1.1 for type-1 AGNs and a = 1.5 for type-2. The steeper slope for 
type-2 was chosen on account of the fact that these dominate over type-1 at low luminosities. As in the 
case of normal late-type and of starburst galaxies, the other parameters are obtained by minimum \ 2 fits, as 
detailed in Sect. [51 and the best fit values are listed, with their uncertainties, in Table [TJ For type-2 AGNs 
pure density evolution was found to be sufficient to account for the data. 



For type-1 AGNs we adopted the mean QSO SED by Richards et al. ( 20061 ). extended to mm wave- 
lengths as described in Sect. 12.31 while for type-2 AGNs we adopted the SED of the local AGN dominated 



ULIRG Mrk 231, taken from the SWIRE library (jPolletta et al. 1 120071) . These SEDs are shown in Fig. [5] 



where the SED of type-3 AGNs associated to dusty star-forming proto-spheroidal galaxies is also plotted 
for comparison. The SED of type-3 AGNs is the most obscured at optical/near- IR wavelengths due to the 
effect of the dense, dusty interstellar medium of the high-z host galaxies. This means that the counts at 
optical/near-IR wavelengths are dominated by type-1 AGNs with type-2 AGNs becoming increasingly im- 
portant in the mid-IR. The 3 AGN populations have appro ximately the same rat i o betw een the rest-frame 
12 //m and the bolometric luminosity, as first pointed out by Spinoglio fc Malkan ( 1989h . 



The type-1 /t ype-2 space density ratio yielded by the model increas es with luminos ity, consistent with 
observations (e.g. lBurlon et al.ll201lh and with the receding torus model (jLawrencdl 1 99 lh . In the framework 
of the standard unified model of AGNs type-1 and type-2 AGNs differ only in terms of the angle which the 
observers line of sight makes with the axis of a dusty torus. If the line of sight to the central region is blocked 
by the torus, the AGN is seen as a type-2. According to the receding torus model the opening angle of the 
torus (measured from the torus axis to the equatorial plane) is larger in more luminous objects, implying 
that obscuration is less common in more luminous AGNs. Since our model implies that type-1 AGNs (but 
not type-2's) are evolving in luminosity, they become increasingly dominant with increasing redshift. 



4. Source counts and contributions to the background 

The surface density of sources per unit flux density and redshift interval is 

d 3 N(S l/ ,z) _ &{log L„,,z)dL v . d 2 V 
dS v dzd£l L v i In 10 dS v dzdfl ' 

where v 1 = v(l + z), 



the comoving volume per unit solid angle is 

d 2 V c (l+z) 2 D%(z) 
dzdfl ~ H E(z) 



(8) 



_ {l + z)L v , 



(10) 
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and the luminosity distance Z?l and the angular diameter distance £>a are related, in a flat universe, by 

dz' 



1 + z H 



E(z'Y 



(11) 



The differential number counts, i.e., the number of galaxies with flux density in the interval S u ± dS v /2 at 
an observed frequency v per unit solid angle, are then 



d 2 N 
dS^dQ 



{S v 



dz 



$(logV,z) dL v , d 2 V 



(12) 



L y i In 10 dS v dzdVl 

The integral number counts, i.e., the number of galaxies with flux density S v > SV.inf at frequency v per 
unit solid angle, are given by 

d 2 V 



dN K -> <? W 



dz- 



dzdQ 



$(logi^/, z)d\ogL h 



(13) 



where v 1 = (1 + z)v and L„< ,i n f is the monochromatic luminosity of a source at the redshift z observed to 
have a flux density S v ,inf- Counts (per steradian) dominated by local objects (z <C 1) can be approximated 

as 



-.2.5 



d 2 N 



1 1 f 

— — = / $(bg L v , z ~ 0)L 3 J 2 d log L v 
4tt 4V7T 7 



The redshift distribution, i.e. the surface density of sources with observed flux densities greater than a chosen 
limit S'jyjnf per unit redshift interval, is 



d 2 N 
dzdQ 



(,Z, S v > iS^inf ) 



d 3 N 



dS'dzdfl 



rdS' 



(15) 



The steepness of the (sub-)mm counts of proto-spheroidal galaxies and their substantial re dshifts imply 
that their counts are strongly affected by the m agnification bias due to gravitational lensing ( Blain 19961: 
Perrottaet aTlbooi Eool iNegrello et al. Il2007h : 



(i 3 iV 1 en S ed(5' l /, Z) 



dfi 



d 3 N(S v /n,z) dP(fi\z) 



(16) 



d log S u dzdVL J „ d\ogS v dzdVt d/u, 

where dP/d/i is the amplification distribution that describes the probability for a source at redshift z to be 
amplified by factor /i. Here we have approximated to unity t he factor 1/(m) tha t would have appeared on 
the right-hand side, as appropriate for large-area surveys (see lJain fc Lima 1 1201 lh - 



We have computed dP/dfj, using the SISSA model worked out by lLapi et al. I (|2012l ). The differential 
counts including the effect of lensing can be computed integrating eq. (ITBl over z. The effect of lensing on 
counts of other source populations and on proto-spheroidal counts at shorter wavelengths is small and will 
be neglected in the following. 

Interesting constraints on the halo masses of proto-spheroidal galaxies come from the auto- and cross- 
correlation functions of intensity fluctuations. A key quantity in this respect is the flux function, d 2 S^/dzdfl, 
i.e. the redshift distribution of the cumulative flux density of sources below the detection limit S v ,iim 



d 2 S u 
dzdQ 



d 3 N 



dS^dzdQ 



Si/dSi/ 



The contribution of a source population to the extragalactic background at the frequency v is 

d 2 N{S„) 



S u - 



dS v dQ, 



-dS u 



(17) 



(18) 
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5. Determination of the best fit values of the parameters 



A minimum \ 2 approach for estimating the optimum values of the parameters of the physical model 
for proto-spheroidal galaxies and associated AGNs is un feasible because of the lengthy calculations required. 
Some small adjustments compared to earlier versions ( Granato et al. 2004 ; Lapi et al. 20061: Mao et al. 



20071 ) were made, by trial and error, to improve the agreement with observational estimates of luminosity 
functions at z > 1.5. An outline of the model, including the definition of the relevant parameters, is 
presented in Ap pendixlAl The chosen v a lues are listed in Table [ H Discussions of physically plausible range s 
can be found in iGranato et al I faoo4, Icirasuolo et all |2005h T Eapi et al.l (feood ). Ishankar etal~l (|2006h . 



Cook. Lapi. fc Granato I |2009t ) and lFan et al. I |201oh 



On the contrary, the minimum x 2 approach was applied to late-type/starburst galaxies and to reac- 
tivated AGNs. The \ 2 minimiza tion was performed using the r outine MPFIT^\ exploiting the Levenberg- 
Marquardt least-squares method jMore Ill978l iMarkwardt Il2009h . 



The huge amount of observational data in the frequency range of interest here and the large number of 
parameters coming into play forced us to deal with subsets of parameters at a time using specific data for 
each subset. The parameters of the evolving AGN luminosity functions were obtained using: 



• the 5-band local QSO luminosity function of lHartwick fc Schade I (|l990j ). 

• the (/-band QSO luminosity functions at z = 0.55 and 0.85 of Croom et al. ( 20091 ). 



the z < 0.75, 1.24 fim AGN luminosity functions of lAssef et al. I (|201ll) 



the bright end [\og(L 60 / L Q ) > 12] of the local 60 /jm luminosity function of iTakeuchi et al. I (|2003l ). 



the Spitzer AGN counts at 8 and 24 /im of iTreister et al. I (|2006l ). 



The B- and (/-band luminosity functions were used to constrain the parameters of type-1 AGNs (type-2 
being important only at the low luminosity end) while the 1.24 /jm luminosity functions were regarded as 
made by a combination of type-1 and type-2 AGNs, the latter being dominant at low luminosities. 

As for the evolving luminosity functions of "warm" and "cold" galaxy populations we used the following 
data sets: 



the IRAS 60 /im local luminosity function of ISoifer fc Neugebauer I (|1991I) . 



• the Planck 350, 550, and 850 /im local luminosity functions of iNegrello et al. I (|2012f) 

• the Spitzer MIPS counts at 24, 70, and 160 /im of iBethermin et al. I l|20ich . 



• the Herschel PACS counts at 160 /tm of iBerta et al. I (|201 lh . 



the Herschel SPIRE counts at 250, 350, and 500 /tm (jBethermin et al. Il2012bl ). 



The fits of the counts were made after having subtracted the contributions of proto-spheroidal galaxies, which 
are only important at wavelengths > 160 /im. The best-fit values of the parameters are listed in Table [1] 
where values without errors denote parameters that were kept fixed, as mentioned in Sect. [3] 



1 http:/ /purl.com/net/mpfit 
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In comparing model results with observational data the instrumental spectral responses were taken into 
account. This is especially important in the mid-IR because of the complexity of the SEDs due to PAH 
emission lines. The monochromatic luminosity at the effective frequency i/ e g in the observer's frame is given 
by: 

L{v eS ) = J " T{v')L v , {1+z) dv' I lT(y')dv' (19) 

where T(y) is spectral response function and the integration is carried out over the instrumental band-pass. 
When the model is compared with luminosity function data at frequency i>i (in the source frame) coming 
from different instruments for sources at redshift z we use the response function of the instrument for which 
v e s is closest to fi/(l + z). In the case of source counts we use the response function appropriate for the 
most accurate data. 



6. Results 

6.1. Model versus observed luminosity functions and redshift distributions 

The most direct predictions of the physical model for proto-spheroidal galaxies are the redshift-dependent 
SFRs and accretion rates onto the super-massive black-holes as a function of halo mass. During the dust 
enshrouded evolutionary phase the SFRs can be immediately translated into the IR (8-1000 /im) luminosity 
functions of galaxies. As mentioned above, according to our model, the transition from the dust-obscured 
to the passive evolution phase is almost instantaneous and we neglect the contribution of passive galaxies 
to the IR luminosity functions. In turn, the accretion rates translate into bolometric luminosities of AGNs 
given the mass-to- light conversion efficiency for which we adopt the standard value e = 0.1. The SEDs then 
allow us to compute the galaxy and AGN luminosity functions at any wavelength. 

In contrast, the phenomenological model for late-type/starburst galaxies yields directly the redshift- 
dependent IR luminosity functions and that for reactivated AGNs yields the 12 //m luminosity functions. 
Again these can be translated to any wavelength using the SEDs described in the previous sections. 

In Fig. [6] the model IR luminosity functions are compared with observation-based determinations at 
different redshifts. At z > 1.5 the dominant contributions come from the stellar and AGN components of 
proto-spheroidal galaxies. These contributions fade at lower redshifts and essentially disappear at z < 1. The 
model implies that AGNs a ssociated to proto- spheroidal galaxies are important only at luminosities higher 



than those covered by the lLapi et al. I (|2011l ) luminosity functions which therefore have been converted 
to bolometric luminosity functions using their galaxy SED, i.e. neglecting the AGN contribution, so that 
log(£i R /L ) = log(£ioo/i ) + 0.21 and log(Li R /i Q ) = log(L 250 /L Q ) + 1.24. At z < 1.5 "warm" and 
"cold" star forming galaxies take over, "cold" galaxies being important only at low luminosities. Type-2 
AGNs (long-dashed pink lines) may dominate at the highest IR luminosities while type-1 AGNs (long-dashed 
light-blue lines) are always sub-dominant (in the IR). 

The scale on the top x-axis in Fig. |6] gives the star formation rates corresponding to the IR luminosities 

/Lira , / SFR 



and is therefore meaningful only to the extent that the AGN contribution is negligible. Moreover, the 
normalization constant applies to high-z proto-spheroidal galaxies whose IR luminosity comes almost entirely 
from star-forming regions. For more evolved galaxies older stellar populations can contribute significantly 
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to the dust heating (|da Cunha et al.ll2012l) ; therefore Lir is no longer a direct measure of the star formation 
rate and therefore the upper scale has to be taken as purely indicative. 

Observational determinations of luminosity functions are available in many wave bands and for many 
cosmic epochs. The comparison between the model and the observed g-band (0.467 /im) AGN luminosity 
functions at several redshifts is presented in Fig. [3 while the comparison in the J-band (1.24 fim) is shown in 
Fig. [HI The conversion from monochromatic absolute AB magnitude Af\ ab to the corresponding monochro- 
matic luminosity vL v {X) is given by \og(vL v /[Lq]) — — 0AM\j±b — l°g(^/[A]) + 5.530. The contribution of 
type-2 AGNs at z < 1.5 strongly increases from the g- to the /-band. Apart from the low-luminosity portion 
of the J-band luminosity function, very likely affected by incompleteness, the agreement between the model 
and the data is remarkably good. 

The comparisons between the global (stellar plus AGN components) luminosity functions yielded by the 
model and those observationally determined at several redshifts and wavelengths are shown in Figs. [5] [11] In 
Fig. [12]we compare model and observed redshift distributions at various wavelengths and flux density limits. 
The comparisons for all the other wavelengths for which estimates of the luminosity function are available 
can be found in the Web site |http://people.sissa.it/~zcai/galaxy^agn/| 

Note that a substantial fraction of sources have only photometric redshifts. Fo r example, the fraction o f 



20101). Only few sources at z > 2 have spectroscopic redshifts ([Berta et al. 1 1201 11 1. Note that photometric 



2010l) 



photometric redshifts is 91% for the VVDS-SWIRE survey with S^^m > 0.4 mJy (IRodighiero et al 
67.5% for the GOODS-N and 36% for the GOODS-S samples with 5 M ».m > 0.0 8 mJy jRodighiero et al 



redshift errors tend to moderate the decline of the distributions at high-z. The effect is analogous to the 
Eddington bias on source counts: errors move more objects from the more populated lower z bins to the less 
populated higher z bins than in the other way. Thus the observed distributions may be overestimated at the 
highest redshifts. In addition, optical identifications are not always complete. On the whole, observational 
estimates of luminosity functions and of redshift distributions may be affected by systematic effects difficult 
to quantify and the true uncertainties may be larger than the nominal values. 



6.2. Model versus observed source counts and contributions to the CIB 



Model and observed source counts at wavelengths from 15 /im to 1.38 mm are compared in Fig. 1131 At 
wavelengths > 350 /im, where, in the present framework, proto-spheroidal galaxies are most important, the 
model provides a simple physical explanation of t he steeply rising portion of the counts, that proved to 



be ve ry hard to account for by ot her both physical ( Havward et al. 20121 Niemi et al. 2012 ; Lacev et al 



20101 ) and phenomenological (e.g. iBethermin et al. Il2012at iGruppioni et al. 11201 II ) models. 



In our model the sudden steepening of the (sub-)mm counts is due to the appearance of proto-spheroidal 
galaxies that show up primarily at z ^ 1.5, being mostly in passive evolution at lower redshifts. Their 
counts are extremely steep because, due to the strongly negative K-correction, the sub-mm flux densities 
corresponding to a given luminosity are only weakly dependent on the source redshift. Then, since the 
far-IR luminosity is roughly proportional to the halo mass, the counts reflect the high-z luminosity function 
whose bright end reflects, to some extent, the exponential decline of halo mass f unction at high masses. This 



situation results in a very strong magnification bias due to gravitational lensing (jBlain lll996uPerrotta et al 



20021 l2003[lNegrello et al. Il2007h . The counts of strongly lensed galaxies depend on the redshift distribution 
of the unlensed ones. Thus, the good agreement between the model and the observed counts of strongly 
lensed galaxies (see the 350 fim, 500 fim and 1380 /im panels of Fig. IT3]) indicates that the model passes this 
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test on the redshift distribution. 



Low-z "warm" and "cold" star-forming galaxy populations become increasingly important with decreas- 
ing wavelength. At A > 160 (im proto-spheroidal galaxies yield only a minor contribution to the counts. The 
AGN (mostly type-2) contribution implied by the model is always sub-dominant. We find a maximum con- 
tribution in the mid-IR. At 15 /im it is ~ 8% up to 1 mJy and then rapidly increases up to ~ 20% above 10 
mJy while at 24 /im it is ~ 7-8% up to 0.5 mJy and increases up to ~ 20% above 2 mJy, in fair agreement 



with the observational estimates (jTreister et al. Il2006t iTeplitz et al. 1120111) . 



Another test on the redshift distribution is provided by the estimated counts in different redshift slices 
(Fig. I14p . although we caution that the true uncertainties may be larger than the nominal ones since the 
observational estimates are partly based on photometric redshifts and on stacking. The consistency between 
the model and the data is reasonably good. 

Figure [TBI shows the contributions of the different populations to the cosmic infrared background (CIB). 
The model accounts for the full CIB intensity over the whole wavelength range. Only at A < 10 /im other 
galaxy populations, such as passively evolving galaxies, become important. According to the model, for 
A > 350 /im the main contribution to the CIB comes from proto-spheroidal galaxies and the fraction con- 
tributed by these objects increases with increasing wavelengths. Below A = 350 /im lower z "warm" galaxies 
take over, with "cold" galaxies adding a minor contribution. AGNs are always sub-dominant. The model 
gives a total (type-1 + type-2 + type-3) AGN contribution of 8.6% at 16 /im and of 8.1% at 24/xm For 
comparison, 



Tcpli tz et al. J 201 lh estimate a contribution of ~ 15% at 16 /zm; iTreister et al. 



Ballantvne fc Papovich (j2007t ) find a contribution of 



(200 



6) and 



10% at 24 /im. It must be noted that these ob- 
servational estimates are endowed with substantial uncertainties: on one side they may be too low because 
strongly obscured AGNs may be missed, on the other side they may be too high because a significant fraction 
of the observed emission may come from the host galaxy. 



7. Clustering properties of dusty galaxies and power spectra of the cosmic infrared 

background 



An important test of our physical model for the evolution of dusty proto-spheroidal galaxies is provided 
by their clustering properties that are informative on their halo masses. A specific prediction of our model is 
that proto-spheroidal galaxies are the main contributors to the CIB at (sub-)mm wavelengths with "warm" 
starburst galaxies becoming increasingly important with decreasing wavelength. Since, in our model, proto- 
spheroidal galaxies are much more strongly clustered than starburst galaxies, the variation in the mixture 
with wavelength translates in quantitative predictions on the frequency dependence of the amplitude of the 
CIB power spectra and on the level of correlations among the maps at different frequencies. 



power spectra obtained by 



Viero et al. 



We have updated the analysis bvlXia et al. I (|2012h taking into account the new auto- and cross-frequency 



([2012 ) from Herschel/ SPIRE measurements and the power spectrum 



at 100 /im derived by iPenin et al 

spectrum at 160 /im. However the amplitude of the latter is anomalously large 



(|2012T) . The latter authors actually give also an estimate of the power 

As an example, for the 



wave-number kg = 0.03arcmin we find that the amplitude normalized to the CIB intensity 

51 _ ^klPjke)] 1 ' 2 



(21) 



[eq. (13) of 



Viero et al 



(|2012ft ] is ~ 0.08-0.09 at 100, 250, 350 and 500 /tm but jumps to ~ 0.2 at 160 /im. 
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Since such a jump over a small wavelength range looks odd we decided not to use the 160 fim power spectrum. 



All the relevant details on the formalism used are given bv lXia et al. I (|2012n . Briefly, the power spectrum 
of the galaxy distribution is parameterized as the sum of the 1-halo term, that dominates on small scales 
and depends on the distribution of galaxies within the same halo, and the 2-halo term, that dominates on 
large scales and is related to correlations among different halos. The Halo Occupation Distribution (HOD), 
which is a statistical description of h ow dark matter hal os are populated with galaxies, is modeled using a 
central-satellite formalism (see, e.g., Zheng et al. 20051) . This assumes that the first galaxy to be hosted 
by a halo lies at its center, while any remaining galaxies are classified as satellites and are distributed in 
proportion to the halo mass profile. The mean halo occupation function of satellite galaxies is parameterized 
as: (iV S at) oc (M v ir/M sa t) QBat , where M v i r is the halo mass and the power-law index a sa t is a free parameter. 
The key parameter in the 2-halo term is the minimum halo mass, M v i r ,min, that determines the amplitude 
of the effective bias function b e e(z). 



In the IXia et al. I (|2012l) paper the only free parameters are the minimum halo mass, M m j njPro t OS ph, 
and the power-law index of the mean occupation function of satellites, a sa t iPro tosphj of proto-spheroidal 
galaxies. This is because the contribution of late-type galaxies to the power spectra at A > 250 /im is 
always subdominant and therefore the parameters characterizing their clustering properties were poorly 
constrained. This is no longer true if we add the 100 /im power spectrum, which, however, still provides only 



weak constraints on a S at,iate-ty P c- We theref ore fixed that param eter to 

^sat,lat< 



-type 



= 1. The fits to the 



(|2012l ) give log(M, 



Herschel/SPIKEi power spectra determined bv lViero et al. 

and a sa t, P rotosph = 1-55 ± 0.05 (1 a errors), close to the values found bv lXia et al. 
do not constrain these parameters further but yield log(M m j I ^i a te-type/-M©) 



rfa) = 12.15±0.04 



(|2012l ). The 100 ^m data 
11.0 ± 0.06. The nominal 

errors on each parameter have been computed marginalizing on the other and correspond to Ax 2 — 1. We 
caution that the true uncertainties are likely substantially higher than the nominal values, both because 
the model relies on simplifying assumptions that may make it too rigid and because o f possible systema tics 



42012)] at 



affecting the data. Our value of M m i n)Pro tosph implies an effective halo mass [eq. (17) of lXia et al. 
z ~ 2 of proto-spheroidal galaxies, making up most of the CIB, M e g ~ 4.5 x 10 12 Mp. This value is close to 



the estimated halo mass of the most effective star formers in the universe. iTacconi et al. I (|2008l ) estimated 
their mean comoving density at z ~ 2 to be ~ 2 x 10~ 4 Mpc~ 3 . For the standard ACDM cosmology this 
implies that they are hosted by dark matter halos of ~ 3.5 x 10 12 Mq. 

The best fit model power spectra are plotted in Fig. [16] where the 1- and 2-halo contributions of proto- 
spheroidal and late-type galaxies are also shown. The relative contribution of the latter galaxy population 
increases with decreasing wavelength and becomes dominant at 100 /im. This trend implies a decrease of 
the level of correlations among the maps with increasing separation in wavelength. As il lustrated by Fig. 1171 



the model is in very good a greement with the c ross- wavelength correlations measured bv lViero et al. 



and defined by [eq. (14) of IViero et al. I ([20121 )] 



C A 



P 



AxB 



xB 



kg kg 



(22) 



8. Summary and conclusions 



Studies of galaxy properties as a function of morphological type (e.g. iBernardi et al. 1120101 ) have high- 
lighted a dichotomy between the luminosity-weighted ages of early- and late-type galaxies. The former are 
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mostly older than 8 Gyr while most of Sb or later-type spirals are younger than 7 Gyr, corresponding to 
a formation redshift z < 1—1.5. Building on this datum we have worked out a model whereby the proto- 
spheroidal galaxies, in the process of forming the bulk of their stars, are the dominant population in the IR 
at z £ 1.5 while late-type galaxies dominate at lower redshifts. The model is 'hybrid' in the sense that it 
combines a physical, forward model for spheroidal galaxies and the early evolution of the associated AGNs 
with a phenomenological backward model for late-type galaxies and for the later AGN evolution. 

To describe the cosmologi cal evolution of proto -spheroidal galaxies and of the associated AGNs we 



adopted the physical model by iGranato et al. I ((2004), upgraded working out, for the first time, the epoch- 



dependent luminosity functions of sources as a whole (stellar plus AGN component), taking into account in 
a self-consistent way the variation with galactic age of the global SED. With only minor adjustments of the 
parameters the model accurately reproduces the observed luminosity functions at all redshifts {z 1.5) and 
IR wavelengths at which they have been determined. The model naturally accounts for the observed positive 
evolution of both galaxies and AGNs up to z ~ 2.5 and for the negative evolution at higher redshifts. This is 
the result of the combination of two competing effects. On one side cooling and free-fall timescales shorten 
with increasing redshift because of the increase of the matter density and this yields higher star formation 
rates, i.e. higher galaxy luminosities at given halo mass. The higher gas densities are also responsible for 
a delay of the AGN switch-off time by feedback implying positive luminosity and density evolution of these 
objects. These effects are thwarted by the decrease in the comoving density of massive halos that prevails 
above z ~ 2.5 causing a decline of the bolometric luminosity functions of both galaxies and AGNs. 

The model also provides a simple physical explanation of the steeply rising portion of the (sub-)mm 
counts, that proved to be very hard to account for by other physical and phenomenological models. The sharp 
steepening is due to the sudden appearance of proto-spheroidal galaxies that do not have, in this spectral 
band, an evolutionary connection with nearby galaxies because their descendants are in passive evolution at 
z ;$ 1.5. Their (sub-)mm counts are extremely steep because, due to the strongly negative K-correction, the 
flux densities corresponding to a given luminosity are only weakly dependent on the source redshift. Then, 
since the far-IR luminosity is roughly proportional to the halo mass, the counts reflect, to some extent, the 
exponential decline of halo mass function at high masses. 

The steepness of the counts imply a strong magnification bias due to gravitational lensing. The counts 
of strongly lensed sources depend on the redshift distribution that determines the distribution of len sing 



optical depths. In fact, this model was the only one that correctly predicted (|Negrello et al. 1 120071 ) the 



strongly lensed counts at 500 urn and the correct redshift dis tribution of bright (S500 (im > 100 mJy) sub-mm 



sources (jNegrello et al. Il2010t iGonzalez-Nuevo et al. II2012I) 



The epoch-dependent luminosity function of late- type galaxies has been modeled in terms of two popula- 
tions, "warm" and "cold" galaxies with different SEDs and different evolution properties. Simple truncated 
power law models have been adopted for the evolution of these populations. "Cold" (normal) late-type 
galaxies evolve (weakly) only in luminosity, while "warm" (starburst) galaxies evolve both in luminosity and 
in density. 

Below z — 1.5 the far-IR emission of proto-spheroidal galaxies and the associated AGNs fade out rather 
rapidly. The AGNs, however, can be reactivated e.g. by interactions. This later phase of AGN emission has 
been described by a phenomenological model analogous to that used for late-type galaxies, distinguishing 
between type-1 and type-2 AGNs. 

In this framework, there is a systematic variation with wavelength of the populations dominating the 
counts and the contributions to the extragalactic background intensity. Above 350 fim the main contributors 
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to the CIB are proto-spheroidal galaxies. In this wavelength range late-type galaxies dominate the counts 
only at the brightest (where normal "cold" star-forming galaxies prevail) and at the faintest flux densities 
(where "warm" starburst galaxies outnumber the proto-spheroids). But these galaxies become increasingly 
important with decreasing wavelength. Proto-spheroids are always subdominant below 250 /im. This strong 
variation with wavelength in the composition of IR sources implies specific predictions for the auto- and cross- 
power spectra of the source distribution, that may help discriminating between different models. Essentially 
all the alternative models have all source populations present over the full relevant rcdshift range. This implies 
a high correlation between the CIB intensity fluctuations at different frequencies. On the contrary, the present 
model predicts a high (close to unity) cross-correlation only at the longest wavelengths (> 500 /im). At shorter 
wavelengths the cross correlation progressively weakens and we expect little cross-correlation between CIB 
fluctuations at, say, 100 and 500 fim. No observational determination is available for correlations among these 
wavelengths, but in the Herschel/ SPIRE wavelength range, where cross correlations have been measured, 
the model results are in good agreement with observations. 

According to our model, the AGN contribution to the CIB is always sub-dominant. It is maximal 
in the mid-IR where it reaches 8.6% at 16/itm and 8.1% at 24 /mi. These contributions are close to, but 
somewhat lower than most observation-based estimates which however are complicated by the difficulty of 
separating the AGN emission from that of the host galaxy. The AGN contribution to the counts is also 
always subdominant. We find a maximum contribu tion in the mid-IR where the model gives AGN fractions 



in fair agreement with the observational estimates (|Treister et al. Il2006t iTeplitz et al. 11201 lh . 



We are indebted to Matthieu Bethermin for several useful clarifications on flux calibration and color 
correction issues, to Roberto Assef for having sent his IR luminosity functions of AGNs in tabular form and to 
Aurelie Penin for having provided a tabulation of CIB power spectra at 100 and 160 /im and clarifications on 
error estimates. We also benefited from useful comments from an anonymous referee. Z.Y.C. acknowledges 
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IHEP. The work has been supported in part by ASI/INAF agreement n. 1/072/09/0 and by INAF through 
the PRIN 2009 "New light on the early Universe with sub-mm spectroscopy" . 



A. Self-regulated evolution of high- z proto-spheroidal galaxies 

The gas initially associated to a galactic halo of mass M V j r , with a cosmological mass fraction /b = 
Mgas/Mvir = 0.165 is heated to the virial temperature at the virialization redshift, z v ; r . Its subsequent 
evolution partitions it in three phases: a hot diffuse medium with mass M- ln t infalling and/or cooling toward 
the center; cold gas with mass M co id condensing into stars; low-angular momentum gas with mass M les 
stored in a reservoir around the central super-massive black hole, and eventually viscously accreting onto it. 
In addition, two condensed phases appear and grow, namely, stars with a total mass A/* and the black hole 
with mass M.. As mentioned in Section [2~T1 we restrict ourselves to the ranges 11.3 < log(M v i r /Af Q ) < 13.3 
and z v ir ^ 1.5. 

The evolution of the three gas phases is governed by the following equations: 
Mi„f = -M cond -M^°, 

Mcoid = A/cond- [l-K(t)]M*-M™ ld -M?s°, (Al) 

Aires = Airflow — A/bh> 
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that link the mass infall rate, Mi„f, the variation of the cold gas mass, M co \d, and the variation of the 
reservoir mass, M res , to the condensation rate of the cold gas, Mcond, to the star formation rate M*, to 
the cold gas removal by supernova and AGN feedback, M^ d and respectively, to the fraction of gas 

restituted to the cold component by the evolved stars, TZ(t), to the inflow rate of cold gas into the reservoir 
around the central super- massive black hole, Mi nnow , and to the back hole accretion rate, Mbh- 



The hot gas cools and flows toward the central region at a rate 



Mcond 



^cond 



with M-* n{ = /b-Mvir and 



t-cond 



where the coefficient is 10% smaller than the value used by iFan et al 



d2010h . 



(A2) 



(A3) 



Note that the cooling and 

inflowing gas we are dealing with is the one already present within the halo at virialization. In this respect 
it is useful to keep in mind that the virial radius of halo (i? v ir — 220(Af v i r /10 13 Mq) 1 / 3 [3/(1 + z v ir)] kpc) is 
more than 30 times larger than the size of the luminous galaxy, and that only a minor fraction of the gas 
within the halo condenses into stars. Indeed, we need strong feedback processes, capable of removing most 
of the halo gas, to avoid an over-production of stars. This implies that any gas infalling from outside the 
halo must also be swept out by feedback; it could however become important for the formation of a disc-like 
structure surrounding the prefo rmed spheroid once it enters the passive evolution phase, with little feedback 
( Cook. Lapi. fc Granato 20091) . As mentioned in Sect. 12.11 the additional material (stars, gas, dark matter) 
infalling after the fast collapse phase that creates the potential well, i.e. during the slow-accretion phase, 
mostly produces a growth of the halo outskirts, and has little effect on the inner part where the visible galaxy 
resides. 



The star formation rate is given by 



M. 



'^cold 



where the star formation timescale is t+ ~ t c 



d/s with s 



(A4) 



5. For a IChabrier I (2003) IMF of the form 



4>(m) = m x with x — 1.4 for 0.1 ^ m ^ IMq and x — 2.35 for m > 1M© we find 1Z ~ 0.54 under the 
instantaneous recycling approximation. 

The gas mass loss due to the supernova feedback is 

with 

A^snesn-Esn „ „ / -^SN 



(A5) 



/?SN = 



0.6 



£bind 

10 51 erg J Vl0 12 M, 



8 x lO- 3 /M y V0.05 



-2/3 



1 + Z 



(A6) 



We adopt the following values: number of SNe per unit solar mass of condensed stars Nsn — 1.4 x 10 2 /Mq 
fraction of the released energy used to heat the gas csn = 0.05; kinetic energy relea sed per SN E'sn 



10 51 ergs; halo binding energy E hind ~ 3.2 x 10 14 (M vir /10 12 M Q ) 2 / 3 ([(l + z)/4] cm 2 s" 2 dMo fc Mao Il2004h . 



The infrared luminosity 



-1000 /im) associated to dust enshrouded star formation is 

L*, IR (t) = fc*,iR x 10 43 ( Af M * r ) erg s' 1 , 
V M Q yr ) 



(A7) 
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where the coefficient fc^jR depends on the SED. We adopt fe^jR ~ 3 ( Lapi et al. 2011 ; Kennicutt 19981) . 

The cold gas inflow rate into the reservoir around the super-massive black hole, driven by radiation 
drag, is given by 

finflow ^ ^(1 - e" rRD ) ~ a RD x 10- 3 M,(1 - e- TRD ), (A8) 



with 



trd(0 



'RD 



\ / Mcold \ / Myir \ 

) U0 12 mJ V10 13 M m / 



-2/3 



(A9) 



.10 12 M©, 

For the strength of the radiation drag we adopt ord = 2.5 and set t^ d = 3.0. The model also follows 
the evolution of the cold gas metallicity, Z co \^(t). An approximate solution of the equations governing the 
chemical evolution is (Lapi et al., in preparation) 

Stj tcond 



^cold 



S7 



-£z{t) 



= (s 7 -l)t/t COI1 



1 



, x ^ 1 I", , x mm (Mz)l 



(AlO) 



where 7 = 1—7?. — /3sn, the metallicity of the primordial inf ailing gas is Z mf = 10 5 , and the mass fraction 
of newly formed metals ejected from stars, £z{t) is given by 



A, 



B z In 



min(t, i 8 atu 







tz 



(All) 



with A z = 0.03, B z = 0.02, t z = 20 Myr, and Saturation = 40 Myr for the Chabrier's IMF (Z Q ~ 0.02). 
Equation (| Al 1|) accounts for the fact that, soon after the onset of star formation, the metal yield, mainly 
contributed by stars with large masses (> 20 M©) and short lifetimes (tz < 20 Myr), is a relatively large 
fraction of the initial stellar mass (/ mc tai > 0.12) while, as the star formation proceeds, it progressively lowers 
to /metal ~ 0.06 as the main contribution shifts to stars with intermediate masses ~ 9 — 20 M© and lifetimes 
tz ~ 20 — 40 Myr, and finally saturates t o values /meta l ~ .013 as stars with masses < 9 M© and long 
lifetimes (Saturation > 40 Myr) take over (jBressan et al.l Il998h . The two parameters Az and Bz depends 
mainly on the IMF. 



The accretion rate into the central black hole obeys the equation 

M BH = min(M™ c , A Edd A/Edd), 



(A12) 



where Mgg c is the accretion rate allowed by the viscous dissipation of the angular momentum of the gas in 
the reservoir 



A/f vlsc 



Mr. 



5 x 10 3 



Kir 



500 km s~ 



Mres 
Mm 



3/2 



1 



Mm 

Mres 



1/2 



(A13) 



with K accr ~ 10~ 2 and V^ ir = GM v 2 / r 3 [4wAyi T (z)p m (z) , A v i r being the overdensity of a virialized halo 
at redshift z Y - lT within its virial radius r v ; r . Msdd = Ai./etEdd is the accretion rate corresponding to the 
Eddington luminosity given the mass to light conversion efficiency e (we set e = 0.1 so that the Salpeter time 
e^Edd = 4.5 x 10 7 yr) and AEdd( z ) is the Eddington ratio that we assume to slightly increase with redshift 
for z > 1.5 

A Edd (z)~0.1(z-1.5) 2 + 1.0 (A14) 
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up to a maximum value AEdd.max = 4. The growth rate of the black hole mass is 

M.(t) = (l-e)M BH (A15) 
starting from a seed mass M* ccd — 10 2 M Q . The bolometric AGN luminosity is 



L. = eM^ = 5.67 x 10«(JL) (— erg 8 - (A16) 



A minor fraction of it couples with the interstellar medium of the host galaxy giving rise to an outflow at a 
rate 



^QSO _ n - /f Mnfxold 

with 



M^cdd = M wind M , (A17) 

Mi n { + JW co ld 



^ISM 

M wind = , (A18) 

^bind 

and 

L§g& - 2 x 10 44 £Q so( ]1 ^ T ) 3/2 erg b" 1 . (AID) 

Lqsq is the mechanical AGN luminosity, used to unbind the gas. The coefficient quantifying the strength 
of the QSO feedback is chosen to be eqjso = 3. The ratio of the mechanical to the total AGN luminosity 



1/2 



L^/L.* 3.5 X10- 3 ^(J^S_] ^ (A20) 



is constrained to be in the range 0.006-0.15. 



Examples of the resulting evolution with galactic age of properties of the stellar and of the AGN 
component are shown in Fig. [T5]for three values of the virial mass and z v ir = 3. 

As mentioned in Sect.[5j to improve the fits of the data we have modified, by trial and error, the values of 
some model parameters used in previous papers, still within their plausible ranges (see Table [2]). The impact 
of these parameters on the derived luminosity functions can be more easily understood wi th reference to th e 



time lag between the halo virialization and the peak in black hole accretion rate, Ai pea k (|Lapi et al. 2006). 
The duration of star formation is A^sf ^ Ai pea k (see Fig. [T8f due to the drastic effect of QSO feedback in 
massive halos which dominate the bright end of the luminosity functions. Note that longer Ai pca k (or A^sf) 
imply higher bright tails of the luminosity functions. The final black hole mass increases with increasing 
the coefficient, t^ d , of the optical depth of gas clouds [eq. (|A9[) ] because it implies a higher efficiency of the 
radiation drag driving the gas into the reservoir. There is a degeneracy, to some extent, between Ten and the 



gas metallicity Z co \d, implying that r^ D cannot be tightly constrained (see lGranato et al. 1120041) . The value 
of Aipoak grows substantially in response to a small increase of the radiative efficiency e that yields a slower 
growth of the black hole mass and a weaker QSO feedback. Higher values of the Eddington ratio, AEdd, 
result in lower values of both A£ pea k and of the final black hole mass. A rise of AEdd at high- z is required 
to account for th e observed space density of very luminous QSOs (see the high-z data in Fig. [7] and Fig. [8] 



Lapi et al.1 120061 ) . A higher QSO feedback efficiency (higher eQso) shortens the duration of star formation, 
AisF, but has a minor effect on Ai pca k and on the final black hole mass. Finally, the coefficient relating the 
SFR to the IR luminosity, fc^jR, varies with age mix of stellar populations, chemical composition and IMF. 
Increasing it we shift the luminosity functions towards higher luminosities. 
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Fig. 1. — Bolometric luminosity functions of proto-spheroidal galaxies. The upper left panel shows the 
luminosity functions at z = 1.5 of the stellar (dot-dashed orange line) and of the AGN component luminosity 
(triple-dot-dashed magenta line), as well as the global luminosity function (solid black line). Note that, as 
discussed in Section [2T2l the latter is not the sum of the two components although, in this case, is very 
close to it. The upper right panel illustrates the evolution of the global luminosity function from z = 1.3 
to z = 4.5, while the lower panels show the evolution of each component separately. The decline at low 
luminosities is an artifact due to the adopted lower limit to the proto-spheroid halo masses. The figure 
highlights the different shapes of the stellar and AGN bolometric luminosity function, with the latter having 
a more extended high luminosity tail, while the former sinks down exponentially above ~ 1O 13 L0. The 
evolutionary behaviour of the two components is qualitatively similar and cannot be described as simple 
luminosity or density evolution; down-sizing effects are visible in both cases. On the other hand there are 
also clear differences. 
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Fig. 2. — SEDs of stellar and AGN components of proto-spheroidal galaxies. The solid black line shows the 
adopted SED for the stellar component, obtained modifying that of the z ~ 2 .3 galaxy SMM J2135-0 102, 



also shown fo r comp arison [solid orange line; the photometric data data are from lSwinbank et al. I (|2010h and 



Ivison et al. 1 l|20ich ]. The dotted magenta line represents the SEP adopted for the dust obscured phase of 



the AGN evolution and is taken from the AGN SED library by iGranato fc Danese I (|1994l ). For unobscured 
AGNs we have adopted the mean QSO SED of Richards et al. (2006; solid magenta line). The original 
SMM J2135-0102 SED and the 2 AGN SEDs are normalized to \og(L m /L Q ) = 13.85 while the modified 
SED is normalized to log(LiR/L Q ) = 13.92 to facilitate the comparison with the original SED. Except in 
the rare cases in which the AGN bolometric luminosity is much larger than that of the starburst, the AGN 
contribution is small at (sub-)mm wavelengths, while it is important and may be dominant, in the mid-IR. 
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Fig. 3. — Global SEDs (solid black lines) for two galactic ages (0.3 and 0.48 Gyr) and three host halo masses 
(log(M^/M Q ) = 11.4, 12.2 and 13.2, from left to right), virialized at z v ; r = 3. The dot-dashed orange line 
(overlaid by the solid black line in some panels) and the triple-dot-dashed magenta line show the stellar 
and the AGN component, respectively. The shorter evolution timescale of the AGNs is clearly visible. The 
effect of feedback as a function of halo mass on the SFR is very different from that on accretion onto the 
super-massive black- hole (see text). 
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Fig. 4. — Adopted SEDs for the "warm" (dashed bl ue line) and "cold" (d otted red line) low-z star- forming; 



galaxies. They were generated combining SEDs of iDale fc Helou I (|2002i ) and ISmith et al. I (|2012l ). as de- 



scribed in the text. The solid orange line shows, for comparison, the SED of proto-spheroidal galaxies. The 
3 SEDs are normalized to the same total IR luminosity log(I/rR/Z©) = 1. 
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Fig. 5. — SEDs of low-z type-1 AGNs (solid light-blue line) and type-2 AGNs (solid pink line). The dotted 
magenta line shows, for comparison, the adopted SED of AGNs associated to dusty proto-spheroidal galaxies 
(type-3 AGNs). The SEDs are normalized to the same, arbitrary, bolometric luminosity. 
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Fig. 6. — Comparison between model and observational determinations of the IR (8-1000 /mi) luminosity 
functions at several redshifts. At z > 1.0 we have contributions from proto-spheroidal galaxies (dot-dashed 
orange lines) and from the associated AGNs (both obscured and unobscured; triple-dot-dashed magenta 
lines). The thin solid black lines (that are generally superimposed to the dot-dashed orange lines) are the 
combination of the two components. These contributions fade at lower redshifts and essentially disappear 
at z < 1. At z < 1.5 the dominant contributions come from "warm" (short-dashed blue lines) and "cold" 
(dotted red lines) star forming galaxies. Type-2 AGNs (long-dashed pink lines) or type-3 AGNs associated to 
dusty proto-spheroids (triple-dot-dashed magenta lines) dominate at the highest IR luminosities while type-1 
AGNs (long-dashed light-blue lines) are always sub-dominant (in the IR). The thick solid black lines show 
the sum of all contributions. The upper horizontal scale gives an estimate of the SF Rs corresponding to IR 



luminosities. These e stimates are on l y indi cative (see Sec t.l6.1|l . Data points are from Le Floc'h et al. ( 20051 



black open squares), Caputi et al. ( 2007 , bl ack stars), Magnelli et al. ( 20091 green do wnward triangles 
Rodighiero et al. ( 2010 . blue open asterisks), Magnelli et al. ( 20 111 black triangles), and Lapi et al. ( 2011 



black open circles). 
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Fig. 7. — Comparison between model and observed g-band (0.467 /j,m) AGN luminosity function at several 
redshifts. As in Fig.[6]the long-dashed light-blue and pink lines refer to type-1 and type-2 AGNs, respectively, 
while the triple-dot-dashed magenta lines refer to AGNs associated with proto-spheroidal galaxies. At 
z < 2 the solid black line shows the sum of all the contributions. At higher z on ly proto-sphero i ds are 



considered. D ata points are from Hartwick fc Schade I ( 1990 , black open circles), Warren et al. 



black crosses), Croom et al 



black open squares) , Palanque-Delabrouille et al. (hold b lack downward triangles) , and lRoss et al 



( 2004, blue stars) . IRichards et al. I ([20061 red triangles') . ICroom et al. 



1994, 



2009 



2012, 



black diamonds). The data bv lrlartwick fc Schade ( 199Clh . given in terms of Mb in the Vega system, were 
converted to M g adopting the B — g ~ 0.14 colour estimated by Fukugit a et al. I (119961) and were further 



corrected for the the different cosmology. The UV magnitudes of lWarren et al. ( 1994) we re first converted 



to B magnitudes {Mb = Mcu3ifc<U ~ t l-39 a K + 0.09, with a v = -0.5) following |PejJ dl99i) and then to M, 



as bef ore. The data bv iRoss et al. I (|2012l ) were converted from Mj,(z — 2) to M g following IRichards et al. 



(2006) with spectral index a v — —0.5. The correction for the different cosmology was also applied. Finally, 
the conversion from M g to vL u (0.467 /jm) is given in Section RTT1 
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Fig. 8. — Comparison between model and observed J-band (1.24 /j,m) AGN luminosity function at several 
redshifts. The lines have the same mea ning as in Fig. [7j There are clear signs of substantial incom pleteness at 



the lowest l uminosities. Data a re from iRichards et al. I (|2006l . red triangles) J Asset et al. I (|201ll black filled 



circles), and lRoss et al. I (|2012l black diamonds). The data by lAssef et al. I ([201 if ), given in terms of Mj in 
th e Vega system were converted to vL v (1.2 41 /urn) using the relatio n (L „(1.241 Mm) = 1623 x io~ - 4Mj Jy) 



by iRieke et al. I ([20081 ) . The z-band data by IRichards et al. I (J2006)) and iRoss et al. I ([20121 ) were converted 



to Mj.ab assuming spectral index a v = —0.5. 
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Fig. 9. — Comparison between mode l and observe d 15 (im global (galaxies plus AG Ns) luminosity function 



at several redshifts. Data are from iPozzi et al. I (|2004l magent a open asterisk s), 



Le Floc'h et al. 



(2005 



black open squares), Matute et al. ( 20061 black filled squ a res), Mazzei et al. ( 2007, black ope n circles), 



Magnelli et al. ( 20091 green triangles). Rodighiero et al. ( 201C, b l ue st ars), Fu et al. ( 2010 . red 



circles for star-f o rmati on and red filled circles for AGNs) , IWu et al. I (|201ll black downward triangles) , and 



open 



Magnelli et al. I (|201ll black triangles). The black filled squares in the panels at z = 0.05 and 0.35 show 
observational estimates of the luminosity function of type- 2 AGNs only while the red filled circles at z = 0.7 
refer to AGN of both types and at z = 1.2 refer to type-1 only. The lines have the same definition as in 
Fig.! 
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Fig. 10. — Comparison between model and observed 90 /mi global (galaxies plus AGNs) luminosity func t ion at 



several redshifts. T he lines have the s a me de finition as i n Fig. [HI Data ar e from ISoifer fc Neugebauer (1991 



asteris ks, 100 /itm), Gruppioni et al. ( 2010l triangles), Sedgwick et al. ( 2011 , diamonds), and Lapi et al. 
(|201lL open circles, 100 /zm) 
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Fig. 11. — Local luminosity functions at (sub-)mm wavelengths. As in the other figures the short-dashed 
blue lines refer to "warm" galaxies, the dotted red lines to "cold" galaxies, the long -dashed pink l i nes to 



type-2 AGNs and the long-dashed l i ght-bl ue lines to type-1 A GNs. Data are f rom Dunne et al. ( 2000l 



orange open squares). I Vaccari et al. I (|201(X light-blue stars), and lNegrello et al. I (|2012[ red open circles). 
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Fig. 12. — Comparison between model and observed redshift distributions at several wa velengths and for 



severa l flux density limits. The l ines have the same defini tion as in Fig. [5] Da ta are from Le Floc'h et al 



(|2009L red open squares, 24 urn). iRodighiero et al. I (120101 blue s tars. 24 jum). iBerta et al. I (|201l[ magenta 
open asterisks, 70, 100 , and 160 /xm), iBethermin et al. I (I2012bl red filled circles, 250, 350, and 500 fixn), 



Chapman et al. I (|2005l . red stars, 850 (an), and lYun et al. I (|2012l . blue open asterisks based on the optical 



photo- z and blue filled circles based on millimetric photo-z). Note that a substantial fraction of sources have 
only photometric redshifts and only few z > 2 redshifts are spectroscopic. Photometric redshift errors tend 
to moderate the decline of the distributions at high-z; thus the observed distributions may be overestimated 
at the highest redshifts (see Section 16. ip . The dip at z ~ 1.5 in the observed redshift distribution of 
sources with Ssso/im > 5mJy is due to the 'redshift desert', i.e. to the lack of strong spectral features 
within the observationa l window and the fast decline at z > 2.5 is due to the lack of radio identifications 
(jChapman et al. 1120051 ). The dip around z ~ 1.5 in the redshift distributions yielded by the model signals 



the transition from the phenomenological approach adopted for low-z sources to the physical approach for 
high-z proto-spheroidal galaxies and associated AGNs. Such artificial discontinuity is a weakness of the 
model that needs to be cleared by further work. 
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Fig. 13. — Euclidean normalized differential number counts at wavelengths from 15 fim to 1380 /im. The 
thick solid lines are the sum of contributions from: "cold" late- type galaxies (dotted red lines), "warm" 
(starburst) late-type galaxies (dashed blue lines), type-1 AGNs (long-dashed light-blue lines), type-2 AGNs 
(long-dashed pink lines), stellar component of proto-spheroids (dot-dashed orange lines), AGN component 
of proto-spheroids (triple-dot-dashed magenta lines), strongly lensed (fi > 2) proto-spheroids (solid green 
lines; only significant at A > 250 /im). The thin solid black lines show the counts of unlensed proto-spheroids, 
including both the stellar and the AGN components; at A > 250 fim these counts essentially coincide with the 
counts of the stellar component only. The filled red circles in the 15 /im panel refer to AGNs only. The filled 



blue c ircles and the open red c i rcles in the 24 /im panel refer to AGNs only and come from iTreister et al 



( 2006 ) and from Brown et al. ( 20061 ) . respect ively. The purple filled squares in the 500 /im panel show the 



estimated counts of strongly lensed galaxies (jLapi et al. 2012). The bright counts at 1.38 mm are also 



generally interpreted as due to strongly lensed galaxies (jvieira et al~ 2010 : Greve et al. 2012 ). References 



for all the data points are given in Table [3] The model provides a physical explanation of the sudden 
steepening of the (sub-)mm counts: it is due to the appearance of proto-spheroidal galaxies that show up 
primarily at z £ 1.5, being mostly in passive evolution at lower redshifts. 
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Fig. 14. — Euclidean normal ized differential number counts per redshift slices. Lines have the same meaning 
as in Fig. [T3l Data are from lLe Floc'h et al. ( 2009, red open circles, 24 /urn). iBerta et al. I (|201lL magenta 
open asterisks, 70 and 100 /an), and iBethermin et al. I (|2012bL red filled circles, 250 and 500 /im). 
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Fig. 15. — Contributions of the different populations to the cosmic infrared background. The lines have the 
same meaning as in Fig. IT3l Proto-spheroidal galaxies are the main con t ributors to the CIB ab o ve ~ 500 ixm. 



Data points are from Renault et al. 1 (2001). Stecker & de Jaeer ( 1997). Laeache et al. 


dl999l). Elbaz et al. 


(2002). Miville-Deschenes et al. (20021). Smail et al. ( 2002 ) . Papovich et al. (2004). Dole et al. (2006). 


Marsdenetal. (20091). 


HoDwoodetal. (2010). Greve et al. (2010). Scott et al. ( 


2010). Altieri et al. 


(2010l). and Berta et al. 


(2011). 
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Fig. 1 6. — CIB angular power spectra at far-I R/sub-mm wavelengt hs. The 100 /im data are from lPenin et al 



([20121 ), those at longer wavelengths are from IViero et al. I (]2012[ ). The lines show the contributions of the 
1-halo and 2-halo terms for the two populations considered here: late-type (LT) "warm" plus "cold" galaxies 
and proto-spheroidal (PS) galaxies. The horizonal dotted magenta lines denote the shot noise level. At 
A > 250 /xm the signal is dominated by proto-spheroidal galaxies while late-type galaxies take over at shorter 
wavelengths. 
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Fig. 17. — CIB cro ss-frequency power spectra at sub- mm wavelengths normalized according to eq. (14 ) of 
Viero et al. I ((2012). The solid line is the result from the model. The data are from lViero et al. I ([2012). 
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Fig. 18. — Evolution with galactic age of properties of the stellar and of the AGN component of proto- 
spheroidal galaxies virialized at z v ; r = 3 for three choices of the virial (mostly dark matter) mass: 
log(M V i r /M Q ) = 11.4 (left-hand column), 12.2 (central column) and 13.2 (right-hand column). In the 
first row, the left y-axis scale refers to masses related to the stellar component [infalling hot gas mass (dot- 
dashed line), cold gas mass (dotted line), stellar mass (M+, solid orange line)] while the right-hand scale 
refers to quantities related to the AGN component [reservoir mass (triple-dot-dashed line) and black-hole 
mass (M, , solid magenta line)] . In the second row the left-hand scale refers to the SFR (dotted line) and to 
the BH accretion rate (dashed black line), while the right-hand scale refer to the IR (8-1000 /im) luminosity 
of the stellar (solid orange line) and of the AGN (solid magenta line) component. In the third row the 
left-hand scale refers to the gas metallicity (solid line) while the right-hand scale refers to the optical depth 
of individual gas clouds (dotted line). 
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Table 1: Parameters for low-z AGNs and for "warm" and "cold" galaxy populations. The parameters of the 
AGN luminosity functions refer to 12 ^im {vL v ) while those for galaxies refer to IR (8-1000 fim) luminosities. 
Values without error were kept fixed. 

AGN 1 (12 /tm) AGN 2 (12 /tm) Warm (IR) Cold (IR) 

log($*/Mpc" a ) -5.409 ± 0.098 -4.770 ± 0.122 -2.538 ± 0.051 -1.929 ± 0.112 



\og{Ll/L Q ) 9.561 ± 0.084 10.013 ± 0.093 10.002 ± 0.076 9.825 ± 0.087 

a 1.1 1.5 0.01 1.372 ± 0.121 

a 0.627 ± 0.017 0.568 ± 0.021 0.328 ± 0.014 0.3 

a* 2.014 ± 0.400 4.499 ± 0.317 0.060 ± 0.200 0.0 

a L 2.829 ± 0.297 0.0 3.625 ± 0.097 1.0 

^brcak 1.0 1.0 1.0 1.0 

Cutoff 2.0 2.0 2.0 2.0 
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Table 2: Parameters of the physical model for the evolution of proto-spheroidal galaxies and associated 
AGNs. The values of the first eight parameters used here are somewhat different from those used in previous 
papers, but still well within the plausible ranges listed in column 3 and discussed in the references given in 
the footnotes. 



Parameter Value Plausible range Description 



RD 


3.0 


1 - 10 a 


Normalization of optical depth of gas cloud [eq. (|A9ll] 


e 


0.10 


0.06 - 0.42 


Black hole accretion radiative efficiency [ea. (]A16(1] 


^Edd 


1 - 4 


< a few b 


Redshift dependent Eddington ratio [eq. (IA14[)] 


e QSO 


3.0 


1 - 10 a 


Strength of QSO feedback [eq. (|Al9|l] 


fc*,IR 


3.1 


2 - 4 C 


Conversion factor from of SFR to IR luminosity [eq. (IA7|)] 


0"* 


0.10 


£ 0.5 


Dispersion of mean stellar luminosity [eq. ©] 


cr. 


0.35 


< 0.5 b 


Dispersion of mean AGN luminosity [eq. ©] 


/gaSjCrit 


0.03 


< 0.165 


Gas mass fraction at transition 

from obscured to unobscured AGNs [§I2.2D] 


esN 


0.05 


0.01 - 0.1 d 


Strength of SN feedback [eq. (jM])] 




2.5 


1 - 10 e 


Strength of radiation drag [eq. (IA8I)] 



Granato et al. 



preparation). 



(2004 


); b 


Lapi et al. 


(2006); c 


Lapi et al. ( 


2011): d 
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Table 3: References for data on number counts (see Fig. [T5|) 



Wavelength (^m) 


Instrument 


Field 


Reference 


15, 24 


AKARI/IRC 


NEP-deep 


Takagi et al. 


(20121 


15 


AKARI/IRC 


de-lensed Abell 2218 




(20101 


15 


AKARI/IRC 


NEP-deep+wide 


Pearson et al. 


(20101 


15 


AKARI/IRC 


CDFS 


Burgarella et al. 


(20091 


15 


AKARI/IRC 


NEP-deep 


Wada et al. 


(20071 


15 


ISO/ISOCAM 


ELAIS-S 


Gruppioni et al. 


C20021 


15 


ISO/ISOCAM 


ISO CAM deep surveys 


Elbaz et al. 


(19991 


16 


Spitzer/IRS 


GOODS-N+S 


Teolitz et al. 


(20111 


24, 70 


Spitzer/MIPS 


ADF-S 


Clements et al. 


(20111 


24, 70 


Spitzer/MIPS 


Spitzer legacy fields 


Bethermin et al. 


(20101 


24 


Spitzer/MIPS 


SWIRE fields 


ShuDe et al. 


(20081 


24 


Spitzer/MIPS 


NDWFS Bootes 


Brown et al. 


(2006) 


24 


Spitzer/MIPS 


GOODS-N 


Treister et al. 


(2006) 


24 


Spitzer/MIPS 


Deep Spitzer fields 


Papovich et al. 


( 20041 


70, 100 


Herschel/PACS 


GOODS, LH, COSMOS 


Berta et al. 


(20111 


70 


Spitzer/MIPS 


xFLS 


Fraver et al. 


(20061 


70 


Spitzer/MIPS 


Bootes, Marano, CDF-S 


Dole et al. 


(20041 


100 


Herschel/PACS 


Abell 2218 


Altieri et al. 


(20101 


250, 500 


Herschel/SPiRE 


HerMES 


Bethermin et al. 


2012b) 


250, 500 


Herschel/SPIRE 


H- ATLAS 


Clement^^^^ 


(20101 


250, 500 


Herschel/SPIRE 


HerMES 


Oliver et al. 


(20101 


250, 500 


Herschel/SPIRE 


HerMES P(D) 


Glenn et al. 


(20101 


250, 500 


BLAST 


BGS P(D) 


Patanchon et al. 


(20091 


500 


Herschel/SPIRE 


H-ATLAS 


LaDi et al. 


(20121 


550, 850 


Planck 


Planck all-sky survey 


Planck Collaboration 


(20121 


850 


SCUBA 


Clusters 


Noble et al. 


(20121 


850 


SCUBA 


Abell 370 


Chen et al. 


(20111 


850 


SCUBA 


Clusters 


Zemcov et al. 


(20101 


850 


SCUBA 


Clusters & NTT-DF 


Knudsen et al. 


(20081 


850 


SCUBA 


SHADES 


Coopin et al. 


(20061 


850 


SCUBA 


Clusters 


Small et al. 


(20021 


870 


APEX/LABOCA 


Clusters 


Johanssoii^^^ 


(20111 


1100 


ASTE/AzTEC 


AzTEC blank-field survey 


Scott et al. 


(20121 


1100 


ASTE/AzTEC 


COSMOS 


Aretxaga et al. 


(20111 


1100 


ASTE/AzTEC 


ADF-S, SXDF & SSA22 


Hatsukade et al. 


(20111 


1100 


ASTE/AzTEC 


GOODS-S 


Scott et al. 


(20101 


1100 


JCMT/AzTEC 


SHADES 


Atistermann et al. 


(20101 


1100 


JCMT/AzTEC 


COSMOS 


Attstermann et al. 


(20091 


1400 


SPT 


SPT survey 


Vieira et al. 


(2010) 



